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Summary 


A  nonlinear  finite  element  analysis  procedure  already  developed  for 
planar  mechanisms  is  being  modified  to  handle  complex  mechanisms  with  sliding 
masses  and  mechanisms  operating  at  relatively  high  speeds.  Progress  is  also 
being  made  in  developing  a  suitable  nonlinear  finite  element  analysis 
procedure  for  three-dimensional  mechanisms.  In  both  cases  the  analysis  takes 
into  account  the  effects  of  geometric  and  material  nonlinearities,  vibrational 
effects  and  coupling  of  deformations.  In  the  optimal  design  area,  a  new 
algorithm  has  been  developed  for  finding  the  minimum  of  a  sum-of-squares 
objective  function  subject  to  general  nonlinear  constraints.  The  solution  of 
preliminary  examples  indicate  good  results  in  terms  of  the  total  number  of 
objective  function  evaluations  to  obtain  an  optimal  design.  Complete  details 
of  these  investigations  are  included  in  the  Appendix.  To  meet  the 
extraordinary  computational  needs  of  this  project,  a  separate  research  VAX 
11/785  Computer  and  peripheral  equipment  were  made  available  through  a  DoD 
research  grant.  The  National  Science  Foundation  also  provided  funds  for  some 
additional  equipment  as  well  as  computational  time  on  a  supercomputer. 


Research  Objectives 


The  objective  of  this  research  is  to  develop  a  nonlinear  finite  element 
dynamic  analysis  procedure  for  planar  as  well  as  spatial  mechanisms  that  are 
frequently  used  in  space  structures.  Included  in  the  nonlinear  analysis  are 
the  effects  of  curvature-displacement  nonlinearity,  nonlinearity  due  to 
extension  or  stretching  (both  caused  by  large  deformations^,  material 
nonlinearity  as  well  as  combinations  of  these.  In  addition  to  the  nonlinear 
analysis,  an  efficient  optimal  design  method  is  to  be  developed  to  handle 
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objective  functions  composed  of  combinations  of  rigid  body  and  deformation 
displacements  involving  geometric  design  variables  as  well  as  cross-sectional 
sizes  of  the  members  of  the  mechanisms  subject  to  limitations  on  stresses  and 
deformations.  During  the  present  reporting  period  it  was  proposed  to: 

1)  Extend  the  nonlinear  analysis  procedure  developed  for  two-dimensional 
mechanisms  to  three-dimensional  mechanisms 

2)  Apply  the  optimization  technique  developed  to  simple  mechanism  problems 

Significant  Accomplishments 

A  major  portion  of  the  research  during  the  reporting  period  was  directed 
towards  achieving  these  research  objectives.  Considerable  progress  has  been 
made  in  both  of  these  areas.  A  summary  of  significant  accomplishments  is 
presented  below  with  complete  details  of  the  work  included  in  the  Appendix. 

A  non-linear  finite  element  analysis  procedure  has  already  been  developed 
for  planar  mechanisms  to  handle  geometric  and  material  nonlinearities. 
Complete  details  of  the  work  are  available  in  a  paper  entitled,  "Finite 
Element  Nonlinear  Vibrational  Analysis  of  Planar  Mechanisms"  by  D.  W.  Tennant, 
K.  D.  Willmert  and  M.  Sathyamoorthy  presented  at  the  ASME  Winter  Annual 
Meeting  in  Miami  Beach,  Florida  in  November  1985  and  published  in  the  Special 
Issue  of  Material  Nonlinearity  in  Vibration  Problems.  The  results  of  this 
investigation  clearly  indicate  the  need  to  include  the  nonlinearities  in  the 
dynamic  analysis  of  mechanisms.  There  are,  however,  a  number  of  difficulties 
in  extending  this  approach  to  problems  of  mechanisms  with  sliding  masses  and 
mechanisms  operating  at  relatively  high  speeds.  To  overcome  these 
difficulties,  the  nonlinear  analysis  procedure  is  being  modified  to  handle 
more  complex  planar  mechanisms.  Considerable  progress  has  been  made  in 
modifying  the  already  developed  nonlinear  analysis  procedure  and  the  required 


computer  programs  for  planar  mechanism  analysis.  Numerical  results  for 
complicated  cases  of  planar  mechanisms,  however,  are  still  not  available.  It 
is  expected  that  new  results  will  be  available  for  presentation  at  the 
AFOSR-SES  (Society  of  Engineering  Science)  meeting  to  be  held  at  the  State 
University  of  New  York  at  Buffalo  in  August  1986.  Research  is  also  in 
progress  to  develop  a  nonlinear  finite  element  analysis  procedure  for 
three-dimensional  mechanisms.  The  formulation  of  these  problems  includes  the 
effects  of  material  and  geometric  nonlinearities  due  to  large  deformations, 
and  vibrational  effects  of  members.  A  complete  report  on  this  will  be  sent  to 
AFOSR  as  soon  as  results  are  available  for  seme  example  problems. 

The  mathematical  optimization  techniques  currently  available  for  solving 
optimal  design  problems  all  require  several  iterations  to  obtain  the  best 
design.  Some  methods  involve  a  large  number  of  iterations,  with  each 
iteration  requiring  numerous  analyses.  These  methods  are  adequate  if  the 
design  problem  is  small,  since  computational  times  are  relatively 
insignificant.  However  for  large  design  problems,  or  ones  in  which  a  complex 
nonlinear  analysis  is  required,  it  is  extremely  important  that  the 
optimization  technique  be  very  efficient,  particularly  in  terms  of  the  number 
of  analyses  required.  In  a  paper  entitled,  "The  Gauss  Optimization  Method  For 
Problems  With  General  Nonlinear  Constraints"  by  T.  E.  Potter,  K.  D.  Willmert, 
and  M.  Sathyamoorthy  (presented  at  the  22nd  Annual  Meeting  of  the  Society  of 
Engineering  Science  held  at  the  Pennsylvania  State  University,  University 
Park,  Pennsylvania,  October  1985),  a  new  algorithm  is  developed  for  finding 
the  minimum  of  a  aum-of-squares  objective  function  subject  to  general 
nonlinear  constraints.  The  solution  of  preliminary  examples  indicate  good 
results  in  terms  of  the  total  number  of  objective  function  evaluations 
required  by  the  algorithm  to  obtain  an  optimal  design.  The  optimization 
techniques  developed  in  this  research,  as  extensions  of  the  Gauss  method  to 
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handle  various  types  of  constraints,  reduce  the  number  of  analyses  required  to 
obtain  an  optimal  design  thereby  reducing  the  computational  time 
significantly.  This  method  is  now  being  used  to  solve  additional  example 
problems. 

Because  of  the  very  complex  nonlinear  analysis  required,  which  must  be 
repeated  many  times  during  the  optimization  process,  a  considerable  amount  of 
computer  time  is  needed  for  this  research.  To  meet  these  needs,  a  proposal 
entitled  "Laboratory  for  Graphical  Analysis  of  Nonlinear  Reformations  in 
Dynamic  St ructural -Mechanical  Systems"  was  submitted  to  DoD  under  the  DoD  - 
University  Instrumentation  Program  to  purchase  a  separate  research  computer 
for  this  project.  This  resulted  in  a  grant  (No.  AF0SR-85-0103)  of  $101,567. 
Although  the  original  proposal  called  for  the  purchase  of  a  VAX  11/730.  a  very 
careful  and  thorough  search  for  the  best  computer  (with  the  available  funds) 
resulted  in  the  purchase  of  a  much  larger  and  faster  VAX  11/785.  Digital 
Equipment  Corporation  offered  a  sizable  reduction  in  cost  of  its  VAX  11/785 
computer  under  the  DEC  Educational  Discount  Program.  Because  of  this 
reduction  and  because  of  additional  cost  sharing  by  Clarkson  University's 
School  of  Engineering,  it  was  possible  to  purchase  the  VAX  11/785  at  no 
additional  cost  to  DoD. 
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Computer  Equipment  Purchased  Through  DoD-URIP  Grant  and  Other  Sources 


I.  VAX  11/785  COMPUTER  HARDWARE 

1.  VAX  11/780  Packaged  System  Including:  $102,750 

(A)  VAX  11/780  CPU 

(B)  2-Mbytes  ECC  MOS  {64-K  chip)  Memory 

(C)  H9652  UNIBUS  Expansion  Cabinet  with  BA11-K  and  DD11-DK 


(D)  VAX/VMS  License  and  Warranty 

2.  TU80  9  Track  Streaming  Tape  Drive  with  Cabinet  8.800 

3.  RUA81  456  Mbyte  Fixed  Disk  19.600 

4.  2-Mbytes  of  Additional  Memory  8,100 

5.  FP780  Floating  Point  Accelerator  8,960 

6.  Two  DMF32-LP  Communication  Interfaces  5,250 

7.  780  to  785  Upgrade  Kit  80,000 

8.  25  ft.  RS232  Sync  Cable  95 

9.  Two  300/1200/2400  Baud  Telephone  Modems  1,060 

10.  Installation  N/C 

11.  Insurance  and  Transportation  1,817 

12.  Miscellaneous  —  Installation  of  Power,  Phone 

Lines,  Terminal  Lines,  etc.  1,756 

Total  Computer  Hardware  Cost  $238,188 

II.  COMPUTER  DISPLAY  TERMINAL 


1. 

Tektronix  M4115B  Computer  Display  Terminal 

$19,950 

2. 

Option  N1 :  Warranty-Plus 

1,025 

3. 

Option  2B:  Additional  512  Kbytes  RAM 

4.600 

4. 

Option  23:  Additional  Four  Planes  of  Display  Memory 

6,000 

5. 

Option  09:  4695  Color  Copier  Interface 

500 

6. 

Option  42:  Single  Flexible  Disk 

1,700 
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7.  Display  Stand  750 

8.  Software  Package  1,000 

9.  4695  Color  Graphic  Copier  1,595 

10.  Option  42:  Warranty-Plus  430 

11.  4926  10  Mbyte  Hard  Disk  4.200 

12.  Option  N1 :  Warranty-Plus  210 

13.  Shipping  371 

Total  Display  Terminal  Cost  $42,331 

III.  SOFTWARE  FOR  VAX  11/785  COMPUTER 

1.  VMS  Operating  System  N/C 

2.  FORTRAN  License  5,170 

3.  DECNET  Communication  Software  2,950 

4.  IGL  Graphics  Software  2,677 

5.  PSI  Access  Software  1,850 

Total  Software  Cost  $12,647 

IV.  TOTAL  HARDWARE  &  SOFTWARE  COSTS  $293,166 


The  total  value  of  the  equipment  and  software  is  $293,166.  Discounts  and 
contributions  from  Digital  Equipment  Corporation,  Tektronix,  Clarkson 
University's  College  of  Engineering  and  the  Department  of  Mechanical  and 
Industrial  Engineering  total  $191,599.  Thus,  the  total  cost  of  the  hardware  and 
software  to  DoD  remained  at  $101,567  as  originally  proposed.  It  should  be  noted 
that  the  capabilities  of  the  VAX  11/785  system,  including  the  Tektronix  4115 
graphic  terminal,  are  enormous  compared  to  the  originally  proposed  VAX  11/730. 
The  VAX  11/785  system  is  5  times  faster  than  the  VAX  11/730,  has  4  Mbytes  of 
memory  (compared  to  only  1  Mbyte  of  memory  for  VAX  11/730),  456  Mbyte  of  disk 
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space  (compared  to  121  Mbyte  of  disk  space)  a  total  of  16  terminal  lines,  and  a 
9  track  streaming  tape  drive  (no  tape  drive  was  included  in  the  original  VAX 
11/730  system) . 

Hardware  and  software  were  also  purchased  to  tie  the  VAX  11/785  into 
Clarkson's  campus-wide  computer  network.  The  physical  link  is  through  the 
University's  VAX  11/780,  but  this  is  tied  to  the  other  computers  on  campus, 
which  is  linked  to  other  universities  though  BITNET.  This  tie  in  of  the  VAX 
11/785  allows  the  users  of  this  research  computer  access  to  many  of  the  other 
facilities  of  the  university,  such  as  high  speed  printers,  digital  plotters, 
laser  printers,  etc.  It  also  allows  researchers  with  terminal  connected  to  the 
other  computers  on  campus  to  sign  on  to  the  VAX  11/785  as  though  they  were 
directly  connected. 

The  Tektronix  4115B  computer  display  terminal,  which  is  connected  to  the 
VAX  11/785  computer,  has  recently  been  expanded  to  improve  its  capabilities. 
Both  a  3-dimensional  wire  frame  and  a  shaded  surface  option  have  been  added. 
These  options  allow  the  terminal  to  locally  manipulate  3-dimensional  objects, 
such  as  rotating  them  in  3-dimensional  space,  removing  hidden  lines,  drawing 
shaded  surfaces,  etc.  These  expansions  result  in  this  terminal  being  equivalent 
to  a  Tektronix  4129  terminal,  which  is  the  most  recent  high-end  terminal 
introduced  by  Tektronix.  The  total  cost  of  these  options  was  $16,475,  which  was 
made  possible  through  contributions  from  Gleason  Foundation,  Proctor  and  Gamble, 
Tektronix,  Corning  Glass  Works  as  well  as  the  University's  School  of 
Engineering. 

A  recent  grant  from  the  National  Science  Foundation  (Grant  No. 
DMC-8500627 ) ,  with  M.  Sathyamoorthy  and  K.D.  Willmert  as  principal 
investigators,  has  included  funds  totaling  $10,715  for  the  further  expansion  of 
the  graphic  facilities.  A  Tektronix  4692  color  graphics  copier  has  been 
purchased  and  a  Tektronix  4107  low  resolution  graphic  terminal  is  anticipated  as 
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a  result  of  this  grant.  These  purchases  should  complement  the  high  resolution 
Tektronix  4115B  terminal  obtained  through  the  DoD-University  Research 
Instrumentation  Program.  In  addition  to  these  equipment  funds,  this  NSF  grant 
provides,  as  part  of  its  Cooperative  Program  on  the  use  of  supercomputers, 
twenty-five  hours  of  CPU  time  on  a  Cray  X-MP  supercomputer. 

A  summary  of  funding  sources  including  the  1984-85  DoD-URIP  Grant  to 
purchase  the  VAX  11/785  research  computer  and  other  associated  equipment 
(including  upgrades)  is  given  below: 


DoD  -  University  Research  Instrumentation  Program  $101,567 
Digital  Equipment  Corporation  Contribution  120,852 
Tektronix  Discount  and  Contribution  11,110 
National  Science  Foundation  10,715 
Clarkson  University's  College  of  Engineering  Contribution  56,999 
Department  of  Mechanical  &  Industrial  Engineering  3,176 
Clarkson's  Educational  Resource  Center  3,300 
Gleason  Foundation  7,000 
Proctor  and  Gamble  5,300 
Corning  Glass  Works  33 7 

TOTAL  $320,356 


As  a  result  of  contributions  from  all  of  these  sources,  the  total  value  of 
the  equipment  within  this  laboratory  exceeds  $300,000  for  an  investment  of  only 
slightly  over  $100,000  from  DoD. 
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ABSTRACT 

A  finite  element  approach  ia  presented  in  this  paper  for  the  nonlinear  vi¬ 
brational  analysis  of  planar  mechanisms.  The  analysis  takes  Into  account  the 
effects  of  aaterlal  and  geometric  nonllnearltles  on  the  dynamic  behavior.  The 
geometric  nonllnearltles  included  in  this  study  are  due  to  stretching  of  the 
neutral  axis  and  the  curvaturc'dlsplacement  nonlinearity,  both  caused  by  large 
deformations.  The  material  nonlinearity  is  due  to  a  nonlinear  stress-strain  re¬ 
lationship  of  the  Ramberg-Osgood  type.  The  analysis  presented  here  makes  use  of 
hermite  polynomials  which  ensure  compatibility  of  curvature  between  elements. 
Using  a  variable  correlation  table,  a  global  system  of  nonlinear  equations  are 
derived  in  terms  of  the  global  unknowns  and  the  kinematics  of  the  mechanism  .  A 
harmonic  series  technique  is  then  used  to  obtain  the  steady  state  solutions  to 
this  system  of  nonlinear  equations.  Numerical  results  are  presented  for  an  ex¬ 
ample  mechanism  and  the  effects  of  the  nonllnearltles  are  discussed. 

INTRODUCTION 

The  importance  of  flexibility  of  linkages  on  the  performance  of  high  speed 
minimum-mass  mechanisms  is  well  recognised.  A  considerable  amount  of  research 
has  been  done  in  this  area  in  the  last  two  decades.  While  it  is  desirable  to 
develop  analytical  and  numerical  procedures  that  enable  the  design  of  rigid  link 
mechanisms  and  robots  to  perform  a  given  function  with  specified  reliability,  it 
is  also  important  to  evaluate  the  effects  of  flexibility  of  elastic  members  on 
on  tnelr  performance.  It  is  known  that  a  mechanism  designed  for  operation  at 
low  speeds  may  not  perform  satisfactorily  at  high  speeds  due  to  the  effects  of 
large  inertia  forces  and  resulting  elastic  deformations.  Thus  it  becomes 
necessary  to  Include  in  the  dynamic  analysis  of  mechanisms,  not  only  the  effect 
of  the  rigid  body  motion,  but  also  the  flexibility  of  the  linkages. 

Most  of  the  previous  investigations  in  the  area  of  elastic  analysis  of 
mechanisms  have  been  carried  out  within  the  framework  of  the  linear  theory  [1- 
17).  However,  Vlscoml  and  Ayre  [18]  used  a  Galerkin-type  nonlinear  analysis 
procedure  to  study  the  vibrations  of  a  slider-crank  mechanism.  A  later  work  by 
Sadler  and  Sandor  [19]  used  the  lumped  parameter  approacn  to  a  nonlinear  dynamic 
model  of  an  elastic  linkage.  The  mechanism  analysed  in  this  paper  was  a  general 
four-bar  linkage  and  the  analytical  model  included  the  response  coupling  associ¬ 
ated  with  both  the  transmission  of  forces  at  the  pin  Joints  and  the  dependence 
of  the  undeformed  motion  of  a  link  on  the  elastic  motloi  of  other  links.  A 
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finite  element  analysis,  with  the  aid  of  the  piecewise  linear  method  of  Martin, 
was  used  by  Sevak  and  McLarnan  (20]  to  carry  out  the  nonlinear  analysis  of  a 
mechanism.  Further  nonlinear  work  dealing  with  the  vibrations  of  elastic  mecha¬ 
nisms  are  reported  in  References  l 21 -23 | ,  In  a  recent  investigation,  Thompson 
and  Sung  [24]  used  a  variational  formulation  for  the  nonlinear  finite  element 
analysis  of  planar  mechanisms  considering  geometric  nonlinearities.  Some  experi¬ 
mental  results  were  also  presented. 

This  paper  is  concerned  with  the  nonlinear  vibrational  analysis  of  general 
planar  mechanisms.  A  finite  element  method  is  used  which  includes  the  effects  of 
both  geometric  and  material  nonlinearities.  The  geometric  nonlinearities  in¬ 
cluded  in  this  study  are  due  to  stretching  of  the  neutral  axis  with  partially 
constrained  ends  and  a  general  curvature-displacement  relationship,  both  caused 
by  large  deformations.  The  material  nonlinearity  is  of  the  Ramberg-Osgood  type 
with  three  parameters  to  represent  the  nonlinear  9tress-strain  relationship  (25- 
27].  Additional  effects  considered  are  transverse  shear  and  rotatory  inertia 
and  changes  in  cross-section  due  to  realistically  proportioned  members.  The 
governing  nonlinear  differential  equations  are  derived  for  each  element  in  terms 
of  the  axial  and  transverse  deformations,  rotations,  curvatures,  and  shear  defor¬ 
mation  angles.  These  equations  are  then  assembled  with  the  aid  of  a  variable 
correlation  table  and  the  resulting  global  system  of  equations  is  solved  using  an 
iterative  technique  based  on  a  harmonic  series  solution  procedure. 

FINITE  ELEMENT  FORMULATION 

A  finite  element  method  is  presented  below  for  the  nonlinear  analysis  of  a 
general  closed  looped  mechanism.  The  mechanism  can  be  composed  of  various  com¬ 
binations  of  simple  four  bar  chains,  frame  elements,  sliders  moving  on  fixed  ref¬ 
erences,  or  sliders  moving  on  rotating  links.  Each  link  is  divided  into  one  or 
more  elements  with  each  element  having  the  local  coordinate  system  as  shown  in 
Figure  l.  If  a  slider  is  present,  the  masses  M^  and  M2  are  located  at  ends  1  and 
2,  as  shown.  The  length  of  the  element  A  is  constant  except  for  links  with 
sliders  moving  along  them. 

The  displacement  vector  of  any  point  (a)  on  the  element's  neutral  axis  is 
given  by: 

S  ■  (Xj  cosy  +  Yi  slny  +  x  +  u)i  +  (Y^  cosy  -  siny  +  w)j  (1) 

where  X^  and  Yj  are  the  coordinates  of  end  1  of  the  element  given  by  the  rigid 
body  motion.  The  coordinate  x  is  measured  along  the  element's  neutral  axis  from 
l  to  2  and  y  is  the  angle  between  the  rigid  body  position  and  the  X-axis.  The 
axial  and  transverse  displacements  of  point  (a)  from  the  rigid  body  position  are 
given  by  u  and  w,  respectively.  This  equation  takes  into  account  both  the  rigid 
body  motion  and  the  elastic  displacements  and  defines  the  position  of  any  point 
along  the  neutral  axis. 

Differentiating  Equation  1  with  respect  to  time  yields  the  velocity  of  any 
point  (a).  The  unit  vectors  i  and  j  move  with  coordinate  system  and  vary  with 
time.  The  angular  velocity  of  any  differential  line  segment  on  the  neutral  axis 
of  the  element  Is  given  by: 

Tf.t  ♦  “.*t  (2) 

where  y,t  is  the  derivative  of  y  with  respect  to  time  and  w,xt  is  the  derivative 
of  the  transverse  displacement  with  respect  to  the  local  coordinate  x  and  time  t. 

The  kinetic  energy  due  to  rotation  of  the  element  is  giver  by: 

h  r  1 

K.E.r  -  |  f*  JA[pAx;S,t|2  +  0lz(y,t  ♦  v,xt)2Jdx  dy 

'2 


*  7  +  i  M2!S>t|2.a  (3) 

where  p  Is  the  mass  density  and  Ax  and  Iz  are  the  cross  sectional  area  and  moment 
of  inertia  of  the  element  respectively.  The  tprm  plz(y,t  +  w«xt^^2  represents 
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the  effect  of  rotatory  inertia.  The  kinetic  energy  due  to  beam  bending  associ¬ 
ated  with  transverse  shear  128)  is: 


K.E.B  9  \  J2  !  {p Iza,J}dx  dy 


where  a  is  a  measure  of  the  transverse  shear  angle. 
The  strain  energy  for  the  element  is  given  by: 


(4) 


U  *  f vol  f  a  dE  dvo1  +  J  wol 
o 


1 

2  T  xy  Y  Xy 


dvol 


(5; 


where  o,  f ,  iXy»  an<*  Yxy  are  normal  stress,  normal  strain,  shear  stress  and 
shear  strain,  respectively.  Por  a  nonlinear  material  of  the  Ramberg-Osgood  type 
(25-27),  the  relationship  between  stress  and  strain  is: 


o  *  Ac 


B  em 


(6) 


where  A  corresponds  to  Young's  modulus  E,  and  Bc°  represents  the  nonlinear  term. 
A,  B  and  tn  are  constants  for  the  particular  material  being  considered.  The  above 
relationship.  Equation  6,  is  valid  only  for  positive  strains.  If  the  strain  is 
negative,  the  following  expression  is  used: 

a  =  A  e  +  B  (-e)a  if  e  <  0  (7) 

The  change  in  sign  of  the  nonlinear  term  results  in  the  same  overall  effect  on 
the  stress-strain  relationship  as  for  positive  strain,  i.e.  either  hardening  or 
softening  depending  on  the  values  of  B  and  m. 

Using  the  Ramberg-Osgood  relationship,  the  following  expression  for  the 
strain  energy  is  obtained  for  positive  strain  e  >  0: 

h  h 

«  -  /h  A  c2  -  =T  8  t"*1)’!*  dv  +  /h  I&  J  cxy  yxy  d*  dy  t») 

1  °  V 

where  x  is  the  axial  coordinate,  y  is  the  transverse  coordinate,  and  CXy  is  the 
shear  modulus.  When  c  <  0  the  equation  is: 

h  h 

U  “  ^h  ^  {2  A  c2  *  irtT  ®  <^)«H'l>dx  dy  ♦ 

~2  °  ’2 

The  nonlinear  expression  for  the  curvature  R 
going  large  deformations  is: 

I  .  _ (10) 

R  (l+w,2)3/2 


J  |  cxy  Yxy  d*  dy 
o 

of  a  planar  static  beam  under- 


Thus  the  strain  is 


e  .  -  1 

R 


(l  «■  w. 5)3/2 


(11) 


Combining  the  geometric  nonlinearities  due  to  stretching  of  the  neutral  axis  and 
the  curvature-displacement  nonlinearity,  results  in  the  expressions  for  normal 
and  shear  strain: 
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(12) 


1  .  ?  j.  y 

2  ’  V‘,‘  (1  -  Db  .,2)1/2 


*xy  ■  «.x  +  « 


(13) 


Substituting  these  expressions  for  strain  into  the  strain  energy  Equations 
(8)  and  (9)  produces,  for  £  >  0  or  t  <  Os 

h 

f  2 


S.E.  "  fb  I*  b  (f  *  (».x  +  \  D, 


ito 

2 


(1  +  Db  .,2)3/2 


®k  l  _ 

-  -  B  (+  u*x  +  T  Da  w»x 

m+1  “  X  “  2 

+  2  G*y  ^w,x  +  a)2)dx  <*y 


U  +  D*,  w, 2)3/2 


_)m+l 


(14,15) 


For  equation  (15)  negative  sign  is  applicable  in  all  the  terms  containing  +  signs 
and  la  valid  for  e  <  0.  In  Equations  (12),  (14)  and  (15),  and  D#  are 

tracing  constants  representing  the  effects  of  material  nonlinearity ,  geometric 
nonlinearity  due  to  curvature,  and  geometric  nonlinearity  due  to  stretching  of 
neutral  axis,  respectively.  Each  tracing  constant  is  equated  to  unity  when  that 
particular  nonlinearity  is  being  considered  and  is  equated  to  zero  when  it  is 
not . 

In  order  to  represent  realistically  proportioned  members,  changes  in  cross 
section  are  included.  Each  element  is  divided  into  sections  of  varying  lengths 
with  constant  area.  The  integrations  involved  in  the  element  equations  are 
carried  out  in  a  piecewise  fashion  with  the  area  in  each  section  taken  as  a  con¬ 
stant.  This  procedure  provides  a  reasonable  approximation  of  variable  cross 
sectional  members  without  having  to  resort  to  large  numbers  of  elements. 

The  Eagrangian  function  L  is  defined  sb 
S  »k 

L  '  kii  (*-*‘*  +  k-e'b  ‘  SE)i*  <l6) 

where  is  the  total  number  of  elements  in  the  kth  link  and  S  is  the  number  of 
links  in  the  mechanism.  Substituting  Equations  (3),  (4)  and  (14)  into  Equation 
(16),  the  Lagrangian  L  can  be  expressed  in  terms  of  the  displacements  u  and  w, 
the  shear  angle  a,  and  the  rigid  body  motion. 

Hermite  polynomials  are  used  to  approximate  u,  w,  and  a  in  order  to  satisfy 
the  boundary  conditions  of  various  types  of  mechanisms  easily  and  to  ensure  in¬ 
terelement  compatibility.  The  axial  deformation  u  is  approximated  by  a  linear 
shape  function  given  by 


u(x,t)  -  Uift)  Hj(x)  +  U2(t)  N2(x> 


(17) 


Similarly,  fifth  degree  polynomial  shape  functions  are  used  to  approxisiate  the 
transverse  deformation  V: 


v( x , t )  -  WX(t)  HU(x)  +  0i(t)  H21(x)  4-  ij(t)  H3l(x) 

+  W2(t)  H^2(x)  +  02(t)  H22(x)  +  m2(t)  H32(x)  (18) 

The  shear  angle  a  is  also  approximated  by  a  fifth  degree  polynomial  in  order  to 
make  it  compatible  with  the  transverse  displacement  w.  Therefore  a  is  assumed 
to  bet 

o(x,t)  -  ai (t)  Hjj (x)  +  *i(t)  H21(x)  +  Xi(t)  H3j(x) 

+  a2(t)  H12(x)  +  *2(t)  H22(x)  ♦  X2(t)  H32(x)  (19) 
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where  *  and  \  are  the  first  and  second  derivatives  of  a,  respectively. 

The  Hermite  polynomials  are  given  by: 

N1(x)  •  1  -  e  ;  N2(x)  -  e  (20) 

Hu(x)  -  1  -  10e3  +  15«4  -  6e5  ;  Hl2(x)  -  10e3  -  15e4  +  6*3 

H23(x)  ■  A(e  -  6e3  +  8e4  -  3*3)  ;  H22(x)  ■  A(-4e3  +  7e4  -  3e3) 

H31(x)  -  42(e2  -  3e3  +  3e4  -  e3)/2  ;  H32(x)  -  42(e3  -  2 e4  +  e5)/2  (21) 

where,  e  *  x/A. 

A  transformation  of  coordinates  is  now  introduced  to  change  from  the  moving 
coordinate  system  associated  with  the  elements  to  global  coordinates.  Only  U}t 
V\,  U2  and  W2  need  to  be  transformed.  The  other  coordinates  are  angles  or  de¬ 
rivatives  of  angles  which  are  not  directional  on  the  X,  Y  coordinate  system  used. 
The  transformations  are: 

Ui  ■  Ui  cos*!  -  V\  sin*i  !  U2  ■  U2  cos*2  -  W2  sin*2 

W!  *  U!  sin*!  +  ^1  COfl$l  »  W2  ■  U2  sin*2  +  W2  cos*2  (22) 

For  pin  connections  the  transformation  angles  *!  and  *2  are  set  equal  to  -y  (the 
rigid  body  angle)  which  transforms  the  coordinates  back  to  the  global  coordi¬ 
nates.  For  sliders  moving  on  rotating  links  the  transformation  becomes  more  in¬ 
volved.  In  this  case,  the  deformation  of  the  driver  link  must  be  transformed  to 
correspond  to  the  axial  and  transverse  deformation  of  the  output  link  (14). 

Substituting  the  expressions  from  Equation  (22)  into  Equations  (17)  and 
(18),  the  global  coordinates  for  the  system  are  then: 


q  ■  [Ui  Wj  0!  mi  «!  *!  U2  W2  ®2  “2  a2  *2  *2lT  (23 

The  Lagrangian  function  i«  then  written  in  terms  of  the  transformed  element  co¬ 
ordinates.  Differentiating  the  Lagrangian  with  respect  to  the  element  coordi¬ 
nates,  the  following  element  equations  are  obtained: 

£  <  »  ,  .  .  0  (24 

dt  3q,t 

In  differentiating  the  expressions  for  the  kinetic  and  strain  energies  in  Equa¬ 
tions  (3),  (4)  and  (14),  it  must  be  remembered  that  A,  which  is  the  upper  limit 
of  integration,  is  a  function  of  time.  The  operations  carried  out  in  Equation 
(24)  results  in  a  system  of  nonlinear  element  differential  equations.  Assem¬ 
bling  the  element  matrices  for  the  particular  mechanism  being  solved  results  in 
the  global  system  of  equations: 


The  M,  C,  and  matrices  are  all  functions  of  time.  The  C  matrix  results 
from  the  kinetic  energy  of  the  system.  No  damping  was  included  in  the  formula¬ 
tion  of  the  problem.  The  C  Q,t  term  was  found  to  be  small  and  thus  was  Ignored 
in  the  analysis.  The  matrix  Ke  is  the  linear  portion  of  the  total  stiffness 
matrix.  It  is  a  function  of  rigid  body  motion,  but  Is  not  a  function  of  the 
deformations  Q.  The  matrix  K^,  however,  is  the  nonlinear  portion  of  the  stiff¬ 
ness  matrix.  It  is  a  function  of  the  deformations  Q.  Equation  (25)  is  thus  a 
nonlinear  system  of  differential  equations. 

The  derivation  of  the  finite  element  Equation  (25)  is  based  on  the  assump¬ 
tion  of  positive  strains  c.  If  the  strain  is  negative  e  similar  derivation  is 
possible,  based  on  Equation  (15)  for  the  strain  energy  rather  than  Equation 
(14).  The  only  difference  in  the  resulting  Equation  (25)  is  in  the  stiffness 
matrix.  Wherever  an  e*"*  occurs  in  the  original  formulation,  it  becomes  (-e)*** 
for  negative  strains.  All  other  negative  signs  resulting  from  the  introduction 
of  -c  cancels  out  in  the  differentiations  required.  Thus  in  order  to  handle  both 


positive  and  negative  strains  the  terms  involving  Em"*  in  the  stiffness  matrix 
were  replaced  by  |e|m"3. 

In  order  to  solve  the  nonlinear  system  of  Equation  (25),  an  iterative  ap¬ 
proach  was  used.  First  the  equations  were  solved  using  the  linear  terms  only, 
i.e.  the  Kj,  matrix  was  ignored.  This  was  accomplished  by  setting  all  of  the 
tracing  constants  D^,  D8  and  equal  to  zero.  The  solution  Q  for  the  linear 
equations  was  then  used  to  determine  values  for  the  nonlinear  stiffness  matrix 
Kn-  Equation  (25)  was  then  solved  again  for  new  values  of  Q,  and  the  process  re¬ 
peated.  Experience  showed  that  this  procedure  converged  in  from  3  to  5  itera¬ 
tions.  To  solve  Equation  (25)  for  Q(t),  for  particular  Kn,  a  harmonic  series 
solution  method  was  used  similar  to  that  of  Bahgat  and  Willmert  (14).  This  ap¬ 
proach  overcomes  problems  with  stability,  due  to  the  time  varying  nature  of  the 
matrices,  that  sometimes  result  from  an  eigenvalue  technique.  The  steady  state 
solution  is  obtained  without  adding  artificial  damping.  The  solution,  without 
the  C  Q,t  term,  is  given  by: 

N 

Q(t)  -  l  ( K  -  n2  u2M)-l  cos  nwt  +  Bn  sin  mot)  (26) 

where  u>  is  the  input  crank  speed,  and  ^  and  Bn  are  solutions  to  the  linear  equa¬ 
tions  : 

N  N-l 

FCtfc)  “  cos  nojtfc  +  Bn  sin  nwt^  ;  k  *  0,1,...,2N-1  (27) 

where  B0  is  set  equal  to  zero.  The  values  of  ty  are  the  times  at  2N  equal  time 
increments  per  revolution  of  the  input  crank  given  by: 

tk  -  for  k  ■  0,1 . 2H-1  (28) 

Computational  experience  indicates  that  a  fairly  accurate  solution  is  obtained 
using  only  a  few  terms  in  Equation  (26).  As  the  number  of  terms  increases  the 
components  of  the  matrix  (K  -  n2  <u2K)  grow  and  thus  the  inverse  (K  -  n2  oi2M)"1 
becomes  small.  The  summation  can  therefore  be  truncated  to  reduce  computational 
time . 

The  stress  in  the  links  is  calculated  by  evaluating  the  strains  from  Equa¬ 
tion  (12).  The  stress  can  then  be  determined  at  any  point  in  an  element  using 
Equations  (6)  and  (7).  To  find  the  maximum  stress  in  an  element  the  maximum 
strain  must  be  found.  Setting  the  first  derivative  of  Equation  (12)  to  zero  and 
solving  the  resulting  expression,  the  position  of  the  maximum  strain  is  deter¬ 
mined.  Once  the  location  is  known,  the  maximum  strain  and  stress  can  be  evalu¬ 
ated  . 

The  above  formulation  is  based  on  the  use  of  the  shear  angle  a,  which  is  ap¬ 
propriate  particularly  for  short  members.  For  long  slender  links  this  quantity 
is  not  required.  The  elimination  of  a  reduces  the  size  of  the  problem  consider¬ 
ably  since  the  nodal  deformations  aj,  nij,  A i ,  ^ 2  an<*  *2  wou14  no  longer  be 

present.  For  long  slender  members: 


(29) 


Using  this  expression,  the  equation  for  strain  energy  (14)  for  positive  shear  e 
reduces  to> 


S.E. 
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A  similar  expression  exists  for  negative  strain.  The  kinetic  energy  also  changes 
if  o  is  not  present.  The  energy  associated  with  transverse  shear,  Equation  (4), 
is  eliminated  and  thus  Equation  (3)  represents  the  total  kinetic  energy  of  the 
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element.  Using  a  procedure  similar  to  the  method  outlined  above*  vibrational 
equations  of  the  same  form  as  Equation  (25)  can  be  obtained,  but  they  will  be 
smaller  in  size.  However,  nonlinear  terms  still  exist  in  the  stiffness  matrix 
due  to  the  material  and  geometric  nonlinearities.  The  method  of  solution  is  thus 
identical  to  that  outlined  earlier. 

EXAMPLE  PROBLEMS 

The  following  example  is  presented  to  illustrate  the  method  of  solution. 

The  nonlinearities  due  to  neutral  axis  stretching,  and  complex  curvature-di s- 
placement  and  stress-strain  relationships  are  all  considered.  A  four-bar  link¬ 
age,  as  shown  in  Figure  2,  is  used  as  the  example  with  all  of  the  members  flexi¬ 

ble  and  made  of  the  same  material.  The  data  for  the  mechanism  is;  Length  of 
input  crank  (AB)  -5.0  in.  Length  of  coupler  (BC)  *  11.0  in,  Length  of  rocker 
(CD)  *  10.5  in.  Fixed  distance  (AD)  *  10-0  in.  Cross  section  of  links  a  rec¬ 
tangular,  Height  of  rectangle  *  1.0  in,  and  Width  of  rectangle  ■  0.25  in.  The 
position  of  the  input  crank  is  zero  degrees  at  t  ■  0  and  the  direction  of  rota¬ 
tion  is  counterclockwise.  The  mechanism  is  divided  into  three  elements  with  each 
link  in  the  mechanism  taken  as  an  element.  The  boundary  conditions  are  that  only 
moment  and  shear  terms  exist  for  the  input  crank's  driven  end  (A).  For  the  pin 
connections  between  links,  there  are  deformations,  rotations  and  shear  terms,  and 
for  the  rocker's  fixed  point  there  are  only  rotation  and  shear  terms. 

First,  the  deformations  in  the  mechanism  were  determined  with  the  shear 
angle  a  present.  In  this  case  the  crank  link  was  rotated  at  100  rad/sec.  The 

material  properties,  anproximat ing  aluminum,  were  as  follows:  A  ■  10.87x10^ 

lb/in^,  B  ■  0.8387x10“  lb/in^,  ra  *  3.0,  and  Mass  density  *  0.0002536 
lb-sec^/in^.  Three  separate  procedures  were  used  to  obtain  numerical  results. 
First  the  problem  was  solved  using  the  linear  analysis  method  of  Bahgat  and 
Willmert  1141,  with  E  ■  A.  Next  the  method  of  this  paper  was  used  with  the  trac¬ 
ing  constants  equal  to  zero.  Thus  a  linear  analysis  was  obtained.  Finally  the 
method  was  applied  with  all  tracing  constants  equal  to  one,  i.e.  a  full  nonlinear 
analysis.  A  representative  deformation  Oj  as  a  function  of  crank  position  is 
shorn  in  Figure  3.  This  is  the  horizontal  deformation  of  the  free  end  of  the 
crank  link.  As  can  be  seen,  the  three  curves  are  very  similar.  The  effect  of 
the  shear  angle  a  is  to  increase  the  deformation  slightly.  For  this  slow  speed 
the  linear  and  nonlinear  analyses  were  almost  identical. 

The  same  problem  was  solved  again  at  a  higher  speed  of  200  rad/sec.  The  re¬ 
sulting  deformation  Uj  is  shown  in  Figure  4.  As  can  be  seen,  high  frequency 
oscillations  started  to  appear,  with  greater  separation  between  the  three  analy¬ 
ses.  At  even  higher  speed  these  oscillations  became  more  predominant  to  the 
point  of  instabilities  in  the  motion  at  very  high  speeds. 

The  revised  form  of  the  analysis  equations  was  considered  next,  i.e.  the 
form  without  the  snear  angle  a.  Here  a  crank  speed  of  150  rad/sec  was  used.  A 
comparison  was  made  of  the  effects  of  the  various  nonlinearities  on  the  deforma¬ 
tions  and  stresses  as  compared  to  the  linear  analysis.  Figures  5  and  6  show  a 
comparison  of  the  linear  and  nonlinear  deformations  U2  (the  horizontal  deforma¬ 
tion  of  the  free  end  of  the  output  link)  caused  separately  by  geometric  and 
material  nonlinearities.  Figures  7  and  8  show  the  maximum  stressei  in  the  con¬ 
necting  link  of  the  mechanism.  As  expected,  the  material  nonlineerity  of  the 
Ramberg-Osgood  type  results  in  deformations  which  are  greater  in  magnitude  than 
those  obtained  using  a  linear  elastic  model.  The  maximum  stress  decreased  due  to 
the  presence  of  the  term  Bcm  substracted  from  the  linear  stress  <*xpress  ion . 

The  geometric  nonlinearities  considered,  namely  curvature  displacement  and 
stretching  of  the  neutral  axis,  both  due  to  large  deformations  produced  mixed 
results  with  deformations  reduced  at  some  points  and  increased  at  other  points. 
The  effect  of  the  geometric  nonlinearities  would  be  expected  to  produce  a  stiff¬ 
ening  of  the  members  (28)  of  the  mechanism  and  thus  produce  smaller  deformations. 
The  increased  deformations  in  this  case  might  be  due  to  the  fret  that  the  deform¬ 
ations  are  in  relationship  to  the  entire  mechanism  and  not  Ju.it  to  an  individual 
beam  element. 

CONCLUSIONS 

The  nonlinear  analysis  procedure,  using  a  finite  element  technique,  is  an 
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effective  method  of  calculating  the  ateady  state  deformations  and  stresses  In  a 
mechanism.  Significant  differences  can  occur  between  the  linear  and  nonlinear 
approaches.  This  was  particularly  true  for  the  stresses  in  the  example  consid¬ 
ered  in  this  work.  Research  is  still  needed  on  the  overall  effect  of  the  shear 
angle  a,  and  a  more  complete  picture  of  the  nonlinear  terms  in  the  analysis  would 
be  of  value.  Additional  nonlinear  effect  should  also  be  investigated,  such  as 
the  effect  on  the  translations  of  one  link  due  to  large  deformations  of  Che  other 
links . 
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Fig.  4  Comparison  of  Deformations  U\  at  200  rad/sec 
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ABSTRACT 

A  new  algorithm  is  presented  For  Finding  the  minimum 
oF  a  nonnegative  objective  Function  subject  to  general 
nonlinear  constraints.  This  algorithm,  based  oF  Gauss' 
method  For  unconstrained  problems,  is  developed  as  as 
extension  to  the  Gauss  constrained  technique  For  linear 
constraints.  The  derivation  oF  the  algorithm,  using  a 
Lagrange  multiplier  approach,  is  based  on  the-Kuhn-Tucker 
conditions  so  that  when  the  iteration  process  terminates 
these  conditions  are  automatically  satisFied.  A  Feasible 
design  is  maintained  throughout  the  iteration  process.  The 
solution  oF  preliminary  examples  indicate  excellent  results 
in  terms  oF  the  number  oF  objective  Function  evaluations 
required  by  the  algorithm  to  obtain  an  optimal  design. 


INTRODUCTION 

The  optimal  design  oF  many,  complex  structural  and  mechanical 
systems  is  hindered  by  the  large  computational  times  involved.  Most 
currently  available  optimization  techniques  require  a  large  number  oF 
analyses  to  obtain  the  optimal  design.  For  small  problems,  or  ones 
in  which  the  analysis  is  simple,  these  methods  are  adequate;  however. 
For  large  problems,  or  where  a  time  consuming  analysis  is  required, 
more  eFFicient  optimization  methods  are  needed.  The  goal  oF  this 
research  was  to  develop  such  methods,  particularly  techniques 
applicable  to  mechanical  mechanism  design  where  the  members  are 
deForming  because  oF  high  speed  motion  and  large  external  Forces. 
Computational  times  to  per Form  a  single  analysis  are  enormous  For 
problems  oF  this  type  involving  large  deFormatlons  with  nonlinear 
material  characteristics.  Thus  the  goal  oF  the  methods  developed  was 
to  reduce  the  number  oF  analyses,  even  at  the  expense  oF  Increased 
computational  eFFort  fn  the  optimization  technique  ltselF,  l.e. 
additional  eFFort  in  Finding  new  candidate  design  points  to  analyze. 
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The  methods  developed  take  advantage  of  the  special 
characteristics  of  the  optimizati on  problem,  similar  to  the 
optimality  criterion  techniques.  This  greatly  improves  their 
efficiency.  For  most  mechanism  design  problems,  the  objective 
function  can  be  formulated  as  a  sum  of  squared  quantities  such  as 
the  difference  between  the  desired  performance  and  the  actual 
performance  of  the  mechanism  at  specified  points  during  its  motion. 
Thus  the  techniques  were  developed  specifically  to  handle  problems 
of  this  type,  although  the  methods  are  applicable  to  objective 
functions  which  are  general  sums  of  nonnegative  quantities,  such  as 
weight.  Many  mechanism  problems  have  constraints  which  are  only 
linear  functions  of  the  design  variables.  Thus  a  special  method  was 
developed  for  problems  of  this  type.  Other  problems  have 
constraints  which  are  linear  or  quadratic,  and  another  method  was 
developed  for  this  case.  Some  mechanism  design  problems  have  more 
general  nonlinear  constraints.  Methods  to  handle  these  cases  are 
currently  being  developed. 

All  of  the  techniques  developed  in  this  work  have  been  based 
on  Gauss'  method  C13  which  is  applicable  to  problems  without 
constraints.  Milde  C23  has  shown  this  method  to  be  particularly 
efficient  on  simple  mechanism  design  problems.  The  research 
presented  in  this  paper  has  extended  this  method  to  handle  various 
types  of  constraints  common  to  more  complex  mechanism  design 
problems. 


FORMULATION 

For  an  unconstrained  sum-of— squares  objective  function 

f  <x)  =  ♦  >,  (1) 

where  ♦  is  a  vector  of  linear  or  nonlinear  functions  +  j  thru  ♦  in 
x,  the  Gauss  method ^f or  calculating  the  next  iteration  of  the  P 
design  variables,  xk  +  l,  given  a  current  design,  x^,  is: 

1-1  *  _ 


k+1 


=  xk  -  [j<^k>jT<Jk>]  j<v;<V’ 


(2) 


where 


J<x> 


I  *  I  _ 

7*. <x>  .  .  .  V+_(x) 

l  I1  lP 


-  viT. 


(3) 


It  is  observed  that  only  first  derivatives  sf  the  ♦  functions  are 
required  and  that  the  new  design  point  is  calculated  directly  from 
the  current  design  without  using  a  step  length  determination  with 
associated  one  dimensional  minimization.  Himmelblau  Cl  3  has  shown 
this  method  to  be  very  efficient  for  unconstrained  minimization 
problems. 


This  technique  has  been  extended  to  handle  linear  inequal  ity 
constraints  of  the  form 


gi<x)  -  bjx  -  e.  £  0,  i  -  1, 
as  well  as  equality  constraints 

-a 

h.  <x)  =*  d.x  -  e.  *  O,  i  *  1, 
1  i  x 


,L. 


(4> 
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In  the  derivation  of  the^optimization  method,  the  4^  functions  are 
assumed  to  be  linear  in  x  of  the  form 


(6) 


where  J  is  a  constant  matrix.  However  the  resulting  technique  is 
applicable  to^problems  in  which  the  4^'s  are  general  nonlinear 
functions  of  x. 


At  iteration  k,  the  L  equality  constraints  and  any  of  the 
inequality  constraints  that  are  active  can  be  combined  and  written 
in  the  form 


BTx  -  C  -  O.  (7) 

-k 

If  at  the  next  iteration,  k+1,  the  variables  xfc+1  are  at  the  optimum 
design,  then  the  Kuhn— Tucker  conditions  will  be  satisfied 

Vf(xk+1>  +  BX  =  O  (8) 

BT;k+1  -  C  -  O  (9) 

and 

X  >  O  <io> 


where  X  is  the  vector  of  Lagrange  multipliers.  The  gradient  of  f  is 
given  by 

Vf  (x>  =  2J4(x ) .  <11) 

•a  a 

Expanding  4(x>  in  a  Taylor  series  results  in 

4<xk+j>  “  4<xk>  ♦  j"*" j\<k+^  -  xj  ♦  (higher  order  terms!  (12) 

▲ 

It  is  noted  £hat  the  higher  order  terms  are  equal  to  zero  if  4  is 
linear.  If  4  is  not  linear  then  these  terms  will  be  neglected  and 
the  expansion  is  only  approx! mate^  If  equation  (12)  is  substituted 
into  equation  (11),  evaluated  at  xk+1,  the  result  is 

vfuk+i)  -  4  jTpk+l  -  :k]l.  (i3) 


This  may  be  substituted  into  the  first  Kuhn-Tucker  condition 


equation  <8> ,  and  then  solved  for  x 


k  +  l* 

«k+1  “  *k  "  t2JjTl  [»♦<£,<»  ♦  B*]- 

Plugging  this  equation  for  x.  into  the  second  Kuhn-Tueker 
condition,  equation  (?) ,  yields 

BT£k  -  C  -  BT  [2J JT]  ^2J+  <x ^ )  ♦  Bx]  -  0. 


(14) 


(IS) 


If  the  sang  set  of  constraints  that  are  active  at  x  were  also 
active  at  x^,  then^B  —  C  ■*  0.  Using  this  result,  equation  (15) 
can  be  solved  for  X  asK 


<  16) 


Substituting  this  back  in£o  equation  (14)  and  simplifying  yields  an 
iterative  expression  for  x  .  which  will  give  the  optimum  solution 
if  the  constraints  that  are  active  at  the  optimum  point  (iteration 
k+1)  are  active  at  iteration  ks 


k+1 


xk  ~ 


-  [jjt]  bJbt[jot]  b| 


-1 


[jjt]  \»i(5?k>. 


(17) 


This  expression  is  equivalent  to  that  derived  by  Paradis  and 
Mil Inert  C33  using  a  Gradient  Projection  method  as  the  foundation. 
The  technique  converges  to  the  optimal  design  in  one  iteration  if 
the  objective  function,  f,  is  quadratic  and  the  starting  point  is  on 
the  constraints  which  are  active  at  the  optimal  design.  If  f  is  not 
quadratic,  the  technique  can  still  be  applied,  but  it  will 
generally  require  several  iterations  to  reach  the  optimal  design. 
When  the  technique  terminates,  the  Kuhn— Tucker  conditions  will  be 
satisfied  independent  of  the  form  of  the  objective  function. 

Paradis  and  Willmert  demonstrated  the  efficiency  of  this  method 
by  solving  several  examples.  One  example  presented  was  the  optimal 
design  of  a  four  bar  mechanism  to  generate  a  desired  coupler  point 
path.  The  Gauss  constrained  technique  was  compared  with  the 
Oavidon— Fletcher— Powel 1  method  using  an  interior  penalty  function 
approach  to  handle  constraints.  Using  four  different  starting 
points,  the  Gauss  constrained  method  required  from  23  to  33 
objective  function  evaluations  whereas  the  Davidon-Fletchei — Powell 
method  required  from  209  to  622.  While  not  all  starting  points 
yielded  the  same  optimal  design,  both  method*  reached  the  same  local 
minimum  from  each  starting  point.  Other  examples  also  showed 
considerable  improvement  over  existing  methods. 

The  Gauss  method  has  also  been  extended  to  include  quadratic 
inequality  constraints  or  quadratic  approximations  to  higher  order 
nonlinear  constraints.  In  this  work  the  constraints  are  assumed  to 


have  the  -form 


g^x)  “  ^xA^x  +  B^x  -  £  O,  i  «  1, . . .  ,M.  (18) 

If  at  iteration  k+ 1  there  are  r  active  constraints  (r  £  M> ,  the 
Kuhn-Tucker  conditions  will  be 


+  +  sJs  ' 

—x,  .  ,A  x  +  8  x  —  C  “  0, 
2  k+1  j  k-*-l  j  k+1  j  ’ 


1  *  •  ■  •  »  r 


and 


(19) 

(20) 


X  2  O  (21) 

•there  the  summation  in  equation  (19)  and  the  j  subscript  in  equation 
(20)  ra(tr  to  the  set  of  active  constraints  only. 


Using  a  derivation  similar  to  that  for  linear  constraints, 
substituting  the  expression  for  the  gradient  of  f,  equation  (13), 
^nto  the  first  Kuhn-Tucker  condition,  equation  (19),  and  solving  for 

xk+l  ProducM 


-i 

T  v 

T-*  *  r  A 

2J  J  ♦  >  A  X 

2JJ  x.  -  VfCx.  )  -  >  B  X 

A  J  J1 

k  k  t*  S  4 

J1 

(22) 


This  expression  for  x  in  terms  of  X  is  now  substituted  into  the 
second  Kuhn-Tucker  condition,  equation  (20),  to  obtain: 


r  ■  —  1  r 

2JJT  J  A  X  2JJTx„  - 

A  j  jJ  l  k 


-  £ 

Vf(x  )  -  l  8  X 


k  J  J 

J-1 


’ 

-1 

T  v 

2JJ  «■  >  A  X 

2JJ  x.  -  7f(x.  >  -  >  B  X 

l  A  j 

k  k  j  jj 

2JJT  +  Tax!  |2JJTx.  -  vf(x.  >  -  J  b  x  I 

A J  A  l  k  k  A  J  jJ 


-  ct  -  o. 


i  *  1, . . . ,r. 


(23) 


These  r  nonlinear  equations  1 
and  the  old  design  variables, 
for  the  values  of  X.  thru  X  . 
into  equation  (22)  which  will 


n^terms  of  the  new  unknowns,  Xj  thru  Xr, 
x.  ,  are  solved  by  an  iterative  process 
The  lambda  values  are  then  substituted 
yield  new  values  for  the  design 
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variables.  It  is  observed  that  the  matrix  2JJ  is  the  matrix  of 
second  partial  derivatives,  G,  of  ^he  objective  function  if  it  is 
quadratic.  Thus,  by  replacing  2JJ  in  equation  (22)  and  (23),  this 
technique  becomes  a  modification  of  the  second  order  method  rather 
than  the  Gauss  method. 

At  the  optimum  design  all  constraints  will  either  be  satisfied 
(less  than  zero)  or  active  (equal  to  zero)  and  each  active  constraint 
will  have  a  corresponding  lambda  whose  value  is  greater  than  or  equal 
to  zero.  If,  at  some  iteration  the  set  of  design  variables  yields  a 
violated  constraint,  then  obviously  the  optimum  point  has  not  been 
reached.  In  this  case,  the  newly  violated  constraint  will  be  added 
to  the  set  of  active  constraints  and  the  procedure  allowed  to 
continue.  If  at  some  iteration  the  set  of  design  variables  yields 
all  active  or  satisfied  constraints,  but  one  or  more  of  the  active 
constraints  has  a  corresponding  negative  lambda,  then  the  optimum 
design  has  also  not  been  reached.  The  negative  lambda  implies  that 
the  iteration  process  would  like  to  move  away  from  the  corresponding 
constraint  boundary  toward  the  feasible  region  where  the  constraint 
is  satisfied.  Thus,  the  constraint  is  dropped  from  the  set  of  active 
constraints  and  the  process  allowed  to  continue.  If  more  than  one 
negative  lambda  existed,  then  constraints  are  dropped  one  at  a  time 
starting  with  the  constraint  with  the  most  negative  lambda. 

A  constraint  is  added  to  the  set  of  active  constraints  if  it 
should  become  either  active  (equal _^to  zei;o)  or  violated  (greater  than 
zero)  when  the  step  is  taken  from  xk  to  In  the  case  whei^e  a 

constraint  becomes  violated,  a  line  is  "drawn"  between  xfc  and  «k+1 
and  the  actual  step  is  taken  to  the  farthest  point  along  the  line  so 
that  no  constraints  are  violated.  In  effect,  this  procedure  is  the 
same  as  stepping  back  from  »<k  +  1  toward  x.  until  the  newly  violated 
constraint  is  just  active  (equal  to  zero).  The  constraint  is  then 
added  to  the  set  of  active  constraints  for  the  next  iteration. 

An  example  problem  with  quadratic  constraints  given  by  Boston, 
Willmert  and  Sathyamoorthy  C43  shows  this  method  to  be  very  efficient 
when  compared  to  the  generalized  reduced  gradient  method  (GRG)'.  The 
problem  consisted  of  finding  the  optimal  design  of  a  four-bar 
mechanism  (minimizing  the  coupler  point  path  error  with  respect  to  a 
given  path)  subject  to  several  linear  constraints  on  link  length  and 
movability.  Additionally,  constraints  were  placed  on  the  crank  pin 
to  limit  its  location  to  the  intersection  of  two  circular  (quadratic) 
regions.  The  program  was  run  for  a  four  by  four  matrix  of  problems 
which  included  four  different  starting  points  and  four  different 
conditions  of  the  quadratic  constraints.  For  all  sixteen  runs,  the 
number  of  objective  function  evaluations  for  this  new  method  ranged 
from  8  to  32  (average  was  IS),  while  the  GRG  method  required  from  303 
to  699  (average  was  502)  evaluations. 

The  interesting  information  here  is  that  the  solution  to  this 
quadratical ly  constrained  four-bar  mechanism  used  no  more  objective 
function  evaluations  than  the  linearly  constrained  four-bar  mechanism 
example  by  Paradis  and  Willmert.  While  these  two  examples  are 
necessarily  different,  this  tendency  toward  requiring  similar  numbers 


of  function  evaluations  for  different  classes  of  problems  is  very 
desirable.  The  net  result  is  that  we  now  have  an  optimization 
procedure  for  sum  of  squared  quantities  objective  functions  subject 
to  linear  and  quadratic  constraints  that  not  only  requires  relatively 
few  function  evaluations,  but  seems  to  be  constraint  order 
independent.  Now  the  need  is  to  determine  a  method  which  will  also 
work  for  higher  order  constraints. 

Boston,  et  al .  C43  attempted  to  apply  this  method  to  higher 

order  problems,  but  met  with  mixed  results.  The  problems  encountered 
seemed  to  tie  in  with  the  higher  order  constraints  rather  than  with 
the  higher  order  objective  functions.  There  are  several  limitations 
implicit  in  the  algorithm  which  appear  to  be  the  source  of  the 
problems  encountered.  The  first  limitation  has  to  do  with  the 
application  of  the  constraints,  equation  (18),  to  the  first 
Kuhn-Tucker  condition,  equation  (19).  When  approximating  a  higher 
order  function  by  a  quadratic  Taylor  series  expansion  about  some 
point  x  ,  not  only  is  the  A.  matrix  a  function  of  xQ,  but  so  is  the 
B.  vector  and  the  C.  scalar?  Thus  the  constraint  approximation, 
equation  (18),  shouid  be  written  as 

gi(x)  =  ^xTA^  <xq) x  +  [Bi<x<))j  x  -  C4  (xQ>  <0,  i  =  1,...,M  (24) 

where 

x  =  Xq  +  ax.  .  (25) 

-a  -a  _a 

As  x  approaches  xQ  (or  ax  approaches  0)  then  this  approximation 
approaches  the  exact  value  of  the  constraint.  Thus  as  the  algorithm 
progresses  along  and  constraints  are  added  and  dropped,  the 
constraints  must  be  reapproximated  at  the  latest  design  to  keep  the^ 
step  size  small.  This  can  be  achieved  by  taking  the  new  values  of  x 
as  generated  by  equation  (22)  and  substituting  them  intg  the^ac^ual 
constraint  equations  to  get  improved  values  for  the  A.(x^),  B. (x^), 
and  (xQ)  terms  in  equation  (24)  with  respect  to  the  current  design 
point. 

The  second  limitation  involves  the  second  Kuhn-Tucker  condition, 
equation  (£0),  which  is  used  to  obtain  the  equation  for  the  new 
values  of  X,  equation  (23).  This  is  simply  the  equation  for  the 
active  constraints.  In  the  original  formulation,  an  attempt  was  made 
at  obtaining  a  linear  approximation  in  X  for  this  constraint 
equation.  This  would  allow  equation  (2 3)  to  be  solved  explicitly  for 
X.  However,  failing  this  an  iterative  procedure  was  employed  to  find 
the  values  for  X.  Now  that  an  iterative  process  is  required  there  is 
no  advantage  in  keeping  a  quadratic  approximation  when  the  actual 
constraint  will  work  just  as  well.  Replacing  equation  (20)  with  the 
active  nonlinear  constraint  equations  will  remove  any  errors  due  to 
the  approximation  process. 

The  stepping  back  procedure  for  violated  constraints,  described 
above,  can  also  be  a  source  of  problems.  With  non-convex  programming 
problems  this  procedure  may  lead  to  a  situation  where  the  algorithm 
cannot  move  away  from  a  non-optimum  design.  Because  the  stepping 
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back  procedure  assumed  a  straight  line  path  between  the  two  design 
points,  it  is  possible,  when  backing  out  of  a  newly  violated 
constraint,  to  move  into  the  violated  region  o-f  the  constraint  that 
Mas  active  at  the  beginning  o-f  the  step.  The  procedure  would  then 
step  back  still  further  until  all  constraints  are  satisfied.  It  is 
possible  to  end  up  with  the  same  set  of  active  constraints  as  at  the 

start  of  the  iteration.  In  this  case  the  next  iteration  will  produce 

the  same  design,  which  may  be  non-optimal. 

Two  alternatives  are  readily  apparent  which  may  solve  this 
problem.  The  first  one  is  that  when  a  constraint  becomes  violated, 
repeat  the  step  but  include  the  newly  violated  constraint  in  the  set 

of  active  constraints.  The  second  alternative  i s  to  move  to  the 

point  where  the  constraint  is  violated,  and  then  iterate  from  there 
without  stepping  back.  Of  course,  the  violated  constraint  is 
added  to  the  set  of  active  constraints.  Boston,  et.al.  C43  looked 
into  this  second  alternative  to  some  extent.  They  reported  that  it 
did  not  always  work.  However,  it  is  not  clear  if  it  was  the  "no 
stepping  back"  that  was  the  cause  of  the  problems  or  if  the  second 
order  approximations  to  the  constraints  contributed  to  the  problems. 

In  summary,  the  Gauss  nonlinearly  constrained  technique  is  very 
effective  at  solving  quadratically  constrained  problems.  No  major 
difficulties  appear  to  exist  which  would  preclude  it  from  solving 
problems  with  higher  order  constraints  once  the  modifications 
discussed  above  are  implemented.  This  method  with  the  proposed 
modifications  is  currently  the  leading  candidate  as  the  best  method 
for  solving  highly  nonlinear  mechanism  design  problems. 


RESULTS 

A  verification  of  the  effectiveness  of  the  Gauss  constrained 
method  applied  to  problems  with  quadratic  constraints  is  obtained  by 
solving  the  Rosen-Suzuki  test  problem  C53: 

minimize  F<x>  =  +  X2  +  Zx3  +  x4  ~  1  -  *^:<2  ~  21x3  7x4 


subject  to: 

gjfio  =  x?  +  *2  +  X3  ♦  +  Xj  -  x2  +  x3  -  x4  -  B  <  ° 

g^<x)  =  x2  ♦  2x2  ♦  x2  ♦  2x2  -  x ^  -  x^  -  10  £  0 

gj<x)  ■  2x2  ♦  x2  ♦  x2  ♦  2x  ^  -  5  S  0 

The  optimum  design  for  this  problem  is  at  x  •  [O, l,2,-l] 

Two  versions  of  the  Gauss  nonlinearly  constrained  technique  and 
the  generalized  -educed  gradient  method,  identified  as  GRG.  were  used 


■from  -four  different  starting  points.  One  version  of  the  Gauss 
non linearly  constrained  technique,  identified  as  GNLC,  uses  the 
stepping  back  procedure  and  requires  a  feasible  starting  design  and 
will  always  Maintain  a  feasible  design.  The  other  version, 
identified  an  GNLC. NS,  does  not  use  the  stepping  back  procedure  and 
has  no  requirement  on  the  feasibility  of  the  design  at  any  stage  of 
the  optimization.  The  results  are  summarized  in  Table  I.  It  can 
easily  be  seen  that  the  Gauss  nonlinearly  constrained  technique  is 
much  more  efficient  with  respect  to  number  of  function  evaluations 
than  the  generalized  reduced  gradient  method. 


STARTING^ 

NUMBER  OF 

FUNCTION 

ALGORITHM 

DESIGN,  x0 

ITERATIONS 

EVALUATIONS 

GNLC. NS 

[0,0, 0,0] 

2 

3 

GC.NS 

[1,1, 1,1] 

2 

3 

GC.NS 

[2, 2,  2,  2] 

2 

3 

GNLC. NS 

[0,0, 1/5,0] 

3 

4 

GNLC 

[0,0, /5,0] 

3 

4 

GNLC 

[0,0,  0,0] 

S 

6 

GNLC 

[1,1,1,1] 

5 

6 

GRG 

[0,0, /5,0] 

9 

83 

GRG 

[0,0,  0,0] 

11 

106 

GRG 

[1, 1,1,1] 

11 

133 

GRG 

[2,  2,  2,  2] 

11 

144 

Table  I:  Comparison  of  Algorithms 


CONCLUSIONS 

The  optimization  techniques  developed  in  this  research  as 
extensions  of  the  Gauss  method  to  handle  vai  ious  types  of 
constraints  are  effective  approaches  to  reducing  the  number  of 
analyses  required  to  obtain  an  optimal  design.  As  a  result,  the 
computational  time  for  large  problems  should  be  reduced 
significantly. 
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